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Abstract 

We present a calculation of the strangeness and charmness contents (A^|ss|A^) and {N\cc\N) of 
the nucleon from dynamical lattice QCD with 2 + 1 flavors. The calculation is performed with 
overlap valence quarks on 2+1-flavor domain- wall fermion gauge configurations. The configurations 
are generated by the RBC collaboration on a 24'^ x 64 lattice with sea quark mass ami = 0.005, 
^T ' anis = 0.04, and inverse lattice spacing a~^ = 1.73 GeV. Both actions have chiral symmetry which 

is essential in avoiding contamination due to the operator mixing with other flavors. Nucleon 
propagator and the quark loops are both computed with stochastic grid sources, while low-mode 
substitution and low-mode averaging methods are used respectively which substantially improve 
the signal to noise ratio. We obtain the strangeness inatrix element /^^ = ms{N\ss\N) / Mj^j = 
0.0334(62), and the charmness content /t^ = mc{N\cc\N) /M^ = 0.094(31) which is resolved from 
zero by 3 (T precision for the first time. 



I. INTRODUCTION 

The strangeness and charmness content of the nucleon are of fundamental importance 
to our understanding of the sea quark contribution to nucleon structure. In particular, 
the sea quark contribution to the scalar u/d quark content {N\uu\N) is crucial to further- 
ing our understanding of the pion-nucleon sigma term. In addition to their relevance to 
nucleon structure, the strangeness and charmness content of the nucleon have drawn re- 
cent interest due to their relevance in dark matter searches [l|-l3| . One popular candidate for 
dark matter is a weakly interacting massive particle (WIMP). The scalar (spin- independent) 
effective four-fermion interaction between the WIMP and the quarks, such as the neutralino- 
nucleon scattering in the context of the MSSM, is asq.XXQiQi pJ; l2l- Additionally, the scalar 
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neutralino-nucleus coupling is found to be much larger than the axial- vector coupling [^ . In 
this case, the scalar contribution to the total x-nucleus cross section is 

a=^Zf, + {A-Z)Uf. (1) 

IT 

where rrir is the reduced x-nucleus mass and 

rUN rUg ^-^ ^ VTin 

q=u,d,s ^ Q=c,b,t ^ 

where N = p,n, and 

m, {N\qq\N) 

^^^ - ^^ ^^^ 

mQ{N\QQ\N) 

JTq - I4j 

rriN 

are the contributions for the light (g)and heavy [Q) quarks respectively. From the trace 
anomaly for the nucleon mass 

mQ{l\\QQ\l\)- 

q=u,d,s Q=c,b,t 

and the heavy-quark expansion [Sj], one can relate the glue condensate in the nucleon to the 
heavy-quark condensate. Eq. ([5]) becomes 

27 
mN= Y. rn,{N\qq\N) + -mQ{N\QQ\N) (6) 

q=u,d,s 



mN= Y. ^AN\m\N)+ Y. rnQ{N\QQ\N)-^{N\G,,G,,\N) (5) 



and Eq. ([2]) is written as 
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Q=c,b,t ^ 



(7) 



q=u,d,s 



where fff = 1 —_y\ij=udsfT- This expression is most often used in the analysis for dark 



j-i-uq=u,d,s 
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matter searches [Jj-y]. Since the couphngs a^g and a-sg contain many terms that are pro- 
portional to the quark mass, e.g. through Higgs exchange, we see from Eqs. ([2]) and ([7]) that 
the total spin-independent neutralino-nucleus cross section is mainly proportional to /t, and 
Jt^. Thus, it is important to determine them precisely. 

At energy scales comparable to Aqcd? perturbative calculations of /r^ and fx^ in the 
nucleon are prohibitively difficult due to the nonperturbative nature of quantum chromo- 
dynamics (QCD). However, these sea quark matrix elements are accessible to lattice calcu- 
lations. The first lattice calculations of fr were done with Wilson fermions on quenched 
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lattices and with heavy dynamical Wilson fermion configurations [6H8|. Those calculations 



gave relatively large values for /^^ on the order of ~ 0.19 
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More recent calculations with 






Wilson-Clover dynamical fermions also yield large values [fx^ ~ 0.1-0.46) |9Nll|. whereas 
calculations with fermions incorporating chiral symmetry result in much smaller Jt^ on the 
order of a few percent 
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19|. 



The large value for fj-^ found using Wilson-type fermions is due to the additive renormal- 
ization of the quark mass due to the lattice-spacing-dependent chiral symmetry breaking. 
As a consequence, there is mixing between the uu and dd operators and the ss operator J20| 



12|. This leads to a subtraction term from (A^|ss|A^) which is proportional to the matrix 



element (^N\uu + dd\Ny Since the latter involves the valence contribution, the subtraction 

turns out to be large. For example, it is found that 

^ 2{N\7ss\N) 

{N\uu + dd\N)' ^ ' 

is changed from y = 0.53(12) to y = —0.28(33) after the subtractions for the Nf = 2 
lattice with the non-perturbatively improved clover fermions 20|. Similarly, y = 0.336(3) 
becomes y = 0.059(37) (28) after subtraction for the Nf = 2, 32^ x 64 lattice with Wilson- 
Clover fermions |2]J]. An alternative way of evaluating the strangeness matrix element is to 
apply the Feynman-Hellman theorem and take the derivative with respect to ttt,^ and m'j^ 
instead of the strange quark mass 2j| . This approach avoids the additive mass part of the 
subtraction (l2J and leads to a small fr, = 0.033(16) (4) (2) |22J. 



To avoid the large systematic errors caused by explicit cliiral symmetry breaking with 
Wilson-type fermions, we instead adopt overlap fermions with exact chiral symmetry on the 
lattice for the valence and the quark loop, for which the quark mass receives no additive 
renormaliz at ion . 



In addition to having small 0{a^) discretization errors 23|, |2J], the overlap fermion that 



so be used for the light and charm quarks 



26|. This allows us to calculate both Jt^ 



we use for the valence quarks in the nucleon can a. 
in the loop insertion with small 0{m?a^) error 25 
and /t,. 

For a heavy quark of flavor Q, it is shown |5| that, to leading order in heavy quark 
expansion, the matrix element mQ{N\QQ\N) is related to the glue condensate in the nucleon, 

— o/ 

aQ ^ mQ{N\QQ\N) ^ --1-{N\GG\N). (9) 

If c, 6, and t quarks are all considered to be heavy and the light quark contributions are 
ignored, the combination of the trace anomaly and Eq. ^ predicts 

(yq = —niN = 70 MeV. (10) 

It would be interesting to check this prediction by a direct lattice calculation, since 
the heavy quarks provide a significant contribution to the WIMP-on-nucleon cross section. 
In previous calculations, the charmness content of the nucleon {N\cc\N) = 0.056(27) has 
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been obtained with the MILC HISQ configurations [16|] with a 2 a signal. In this work, we 
calculate the strangeness (A^|ss|A^) and charmness content {N\cc\N) of the nucleon in lattice 
QCD with overlap valence quarks on 2+1-fiavor domain-wall fermion gauge configurations. 

The calculation is done on a 24^ x 64 lattice with sea quark mass mi = 0.005, rris = 
0.04, and lattice spacing a~^ = 1.73 GeV. We obtain the strangeness content fr^ = 
ms{N\ss\N) / M^- = 0.0334(62) which has more than 5 a precision. Compared with pre- 
vious lattice calculations either with disconnected insertion or via the Feynman-Hellman 
theorem, the present calculation has the smallest error. We also obtain the charmness con- 
tent Jt^ = ■mc{N\cc\N) /Mn = 0.094(31). This is the first time the charmness has been 
calculated with more than 3 a precision. 

Our paper is organized as follows: The overlap formulation is briefly summarized in 
Sec. [TTl The technical details for calculating the nucleon two-point functions and the quark 
loops are described in Sec. [1111 and Sec. [IV] respectively. They are combined to calculate 



the disconnected three-point functions in Sec. |Vl Finally, the numerical results with chiral 
extrapolation are presented in Sec. I VI I and the conclusions are given in Sec. IVIII 

II. OVERLAP FERMIONS 

We adopt the overlap fermion formulation for the valence quarks in the nucleon correlation 
functions as well as for the quark loops. The inversion of overlap fermions using deflation of 
low eigenmodes and the construction of meson and nucleon two-point functions with low- 



mode substitution have been detailed previously [26|. Deflation with low eigenmodes and 
HYP smearing speed up the inversion by a factor of ~ 50 for the 24'^ x 64 lattice that we use 
for this calculation and low-mode substitution (LMS) improves the errors of the meson and 



nucleon correlators by a factor of ~ 3 - 4 26| . The overlap quark propagators are calculated 
on gauge configurations with 2 + 1-flavors of dynamical domain- wall fermions (DWF). As we 
mentioned in the introduction, we adopt both the overlap and DWF fermion formalisms since 
they preserve chiral symmetry via the Ginsparg- Wilson relation. Due to the high degree of 
chiral symmetry, the calculation of the quark content is free of the problems that plague the 
Wilson-type fermions as outlined in the above discussion. As an additional advantage, the 
0{im?'a^) discretization errors are small for the overlap fermion. This allows us to compute 
the charmness contribution to the nucleon in addition to the strangeness contribution. 
T^eoveHap ope... 0. defined. 

Do,(p) = l+75£(75/^w(p)), (11) 

where e is the matrix sign function and D„{p) is the usual Wilson fermion operator, except 
with a negative mass parameter —p = 1/2k — 4 in which Kc < i^ < 0.25. We set k = 0.2 in 
our calculation, corresponding to p = 1.5. The massive overlap Dirac operator is defined as 

Tfl Tfi 

= P + y + (P - y) 75 ^(75^«,(p)). (12) 

To accommodate the chiral transformation, it is usually convenient to use the chirally reg- 
ulated field ifj = {1 — ^Dov)^^ in lieu of ip in the interpolation field and the currents. This 
leads to an effective propagator 

G.fl-.(l-%,B-'^^, (13) 
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where Dc = x-d"" 12 ^^ chiral, i.e. {75, Dc} = 28|. It is worthwhile to point out that this 
effective propagator has the same form as that in the continuum 25|. In other words, the 
inverse of the propagator is a chirally invariant massless Dirac operator plus the quark mass 
term. As long as the 0{w?a}^ error is small, this formulation is suitable for both light and 
heavy quarks. 

We adopt the Zolotarev approximation to evaluate the matrix sign function. This entails 
two nested conjugate gradient loops to calculate the propagator of the overlap fermion. For 
each conjugate gradient loop, we use deflation with low eigenmodes to speed up the inversion. 
The details are given in Ref. (26[. Due to the normality of Dqv, i-e. D\^Dov = DovD]^ and 
the Ginsparg-Wilson relation {75,Do,;} = Dovj^Dov, the eigenvalues of Dov are on a unit 
circle with the center at unity. The real and chiral modes are at and 2. Others on the 
circle are paired with conjugate eigenvalues. In other words, if Vi is an eigenvector 



-^ov^i '^i^i: 



(14) 



then its conjugate partner •j^Vi is also an eigenvector with eigenvalue A*, 



Dovl^Vi = A*75fi. 



(15) 



To compute the quark propagator, we first find a few hundred pairs of the lowest eigenvectors 
of the massless overlap operator in addition to the zero modes. Once we have obtained these 
lowest eigenvectors, we can solve the high-mode part of the quark propagator by projecting 
out the low modes from the source 



DeffG 



^V= [^~Y1 [^^^i + 75^*^75) (1 - 2^A„o) J r] 



(16) 



where 77 is the source vector and the factor 1 — I^a^.o takes care of the zero modes which are 
either left-handed or right-handed. 

In contrast, the low-mode part of the effective quark propagator can be constructed with 
eigenvectors directly. 



G- = Y^ 



p\i + m{l-^) pA* + m(l-^,j 



'1 - 2^A,:,o) 



(17) 



For a given source r/, the total effective quark propagator is 



Gr] = G"r] + G^rj 



(18) 



However, we should point out that since G^ is constructed from the eigenmodes rather than 
by inverting a source vector, we can compute the any-to-any propagator for any source and 
sink location with little additional computation. We shall use this fact to carry out the low- 
mode substitution to replace G^r] in Eq. f lTH]) with G^ for a source with given grid points 
which greatly improves the nucleon correlator. This will be explained in more detail in the 
next section. 

III. THE NUCLEON TWO-POINT CORRELATION FUNCTION 

Various attempts have been made to improve the statistics of hadronic two-point correla- 
tion functions such as using a smeared source, a volume source with fixed gauge, and all-to-all 
propagators. The computation of all-to-all propagators usually involves noise sources on dif- 
ferent sites. However, the quark propagator from one site can be contaminated by those 
from neighboring sites. For example, when constructing the nucleon correlation function, 
the three quark propagators may be from the same source site or from different source sites. 
The latter case is not gauge invariant and will introduce noise after averaging over a finite 
number of configurations. Fig. [1] shows the gauge invariant and the non-invariant parts of 
the correlation function. 



The signal-to-noise issue has been examined 26| for the connected hadron correlators 
from the noise source on a time slice and it was found that the noise wall source was worse 
than the point source for all mesons except the pion. It is worse still for the nucleon. In this 
case, the signal to noise ratio is 

cr V V3 

where N is the number of noises and V3 is the 3-volume of the time slice. In addition to the 
usual exponential suppression in time, there is a pre-factor which reduces the signal to noise 
ratio further by a factor of ^/V^. The situation can be ameliorated by reducing the noise 
contamination from neighboring sites with less source points. This introduces the idea of 
a noise grid source with support on some uniformly spaced grid points on a time slice, but 
this does not fundamentally alter the conclusion that, given the same computer time, any 
noise source is worse than the point source for the meson and nucleon. These observations 
suggest a new algorithm for the grid noise with low-mode substitution (LMS) to reduce the 
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FIG. 1. (Color online) Diagram to illustrate the signal and noise of the nucleon correlation 
function with a Z3 noise grid source on a time slice. The upper part with three quarks originating 
from the same spatial site is an example of the gauge invariant signal and the lower one is an 
example of the gauge non-invariant noise with three quarks originating from different spatial sites 
which will be suppressed by gauge average and noise average. 

variance from noise contamination while simultaneously addressing the low-mode correlation 
at the same time 261]. The idea is to replace the low-mode part of the quark propagator 
G^rj estimated from noise sources by the exact G^ in Eq. ( 1T71) . when all three quarks are in 
the low modes or when two quarks are in the low modes and the third is in the high modes, 
in constructing the baryon correlators to reduce the noise contributions shown in Fig. [TJ It 
turns out this LMS is quite successful. It reverses the above-mentioned trend that a noise 
source is worse than a point source and instead reduces the errors compared to the point 



source 



26|. 



The Z3 grid noise is 






(20) 



where ^ is a sparse grid of lattice sites on a time slice, and 6i is the Z^ random phase on 
site i with the property 6f = 1 and (OiOjOk) = 6ij6jk- Tests on the 32^ x 64 lattice with a 
Z3 noise on 64 evenly spaced grid points on a time slice and with LMS reveal that for light 
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quarks (pion masses at 200 - 300 MeV) the errors of the mesons and nucleon masses can be 



reduced by a factor of 3 to 4 as compared to the point source J26|. 

Since the nucleon has a physical size of ~ 0.6-0.8 fm as deduced from its axial and 
electromagnetic form factors, a smeared source for the nucleon usually leads to a reduction 
of errors for the nucleon mass from the point source. In view of this, we introduce a smeared- 
grid source to increase the overlap with the nucleon ground state wavefunction and diminish 
the contribution from the radially excited states and the collateral vrA^ scattering states. We 
adopt the gauge-invariant spatial Gaussian smearing |29| on the grid source. 



r]^{x') = S{x',x)T]{x), (21) 

where S{x',x) is the smearing operator. By design, the smearing operator S{x',x) should 
produce a Gaussian distribution with a Klein-Gordon propagator. It is computed as an 
iteration of many small smearing steps. 



S{x',x) 



f 37/;2\" r 7/,2 JL / 

\ / A — 1 



(22) 



where w is the input width for Gaussian distribution and n is the number of smearing steps. 
The corresponding quark propagator with smeared source is 






B,D-\y,x')S{x\x)b^^,. (23) 



To check if the actual width from the smearing procedure is consistent with the input 
width parameter w, we first define the actual smearing size by 



where, for each spatial position x, p[x) is the norm over spin and color of the smeared source 
vector which is created from a point source vector (5(x). Then we plot r vs w for different 
step sizes n in Fig. [2l We see that for a given ra, there is a range of w where the resultant 
smearing size r has a nearly linear relation with w. Beyond that range, r flattens off as w 
increases. As n increases, the range for the near linear relation expands. Even though we 
don't have evidence that the smeared source, when its r is much smaller than the input ly, 
leads to ill physical effects, we think it is safe to keep the linear relationship. Consequently, 
we will use sufficient n for a prescribed width w so that r and w are in the linear range. 

9 




FIG. 2. The plot shows the relation between the input width parameter w and the output radius 
of the gauge-invariant Gaussian smearing r for different iteration numbers n. 




FIG. 3. The plot shows the profile of the smeared source vector with w = 4: and n = 100 on a 
24^^ X 64 configuration. 
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We choose the parameters w = 4 and n = 100 for this work which give a smearing size 
r ~ 3. The smeared distribution for a typical configuration is plotted in Fig. [3] which is 
indeed close to a Gaussian shape. This particular choice of w and n is aimed at increasing 
the overlap with the nucleon wavefunction while simultaneously minimizing contamination 
from different grid sites. 

In addition to the smearing source, we also placed two sources on the time slices t = and 
t = 32 and calculated the inversion simultaneously. As we shall see later in the calculation 
of nucleon mass and quark scalar matrix elements, the time window for the fitting range 
is between t = 6 to t = 14 which is far from the two sources such that the contamination 
from the time backward propagating 5*11(1/2") state from the second source at a distance 
of 32 time slices away is negligible. This approach nearly doubles our statistics without 
computational overhead. 



The local interpolating operator of the nucleon is taken to be 30 1 



X^,{x) = -e ijy [C)^.p4p, (xjV'^/ [x] (25) 

where C = 727475 in the Pauli-Sakurai gamma matrix convention. The color indices are 
denoted with Latin letters and the Dirac indices are denoted with Greek letters. The nucleon 
correlation function is constructed as 



+trfrG(")'^'^'(2/,x))trfe'')''''(y,x)G(")^'='(i/,x)) ) (26) 



where the G^'^^'^^^'^' {y,x) stands for the u/d quark propagator from the site x to y and the 
color indices from a' to a. The average is taken over the different gauge configurations and 
noise sources. G is defined as {CGC~^)'^. The trace and the transpose operations only act 
on Dirac indices. When the masses of the u and d quarks are set to be equal, G^'^' and G^"^' 
are the same propagator. The correlation function C(G'i, 6*2, Gs) is a linear functional of 
the three functions 6*1,6*2, and G3. 

From Eq. f l23|) . we obtain the quark propagator as a summation of the propagators from 
different grid sites with different Z3 phases: 

G{y,7^') = J2(^^G,{y) (27) 
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where 

G,{y) = D-\y,x')S{x',Xi). (28) 

The nucleon correlation function can be written 

C{G,G,G) = {C{Y,^,G,^Y.^3Gj.Y.^kGu)) 

i j k 

= / ^{.di6j6kG{Gi,Gj,Gk)) 

i,j,k 

= 2_^^ij^jk{C{Gi,Gj,Gk)) 

i,j,k 

^J2iC{G..,G,,G,)), (29) 

i 

where G denotes the correlator in a gauge configuration. G{Gi,Gi,Gi) is the correlator of 
a point or a smeared source at site Xi for a gauge configuration. Thus, given a sufficient 
number of noise vectors and/or gauge configurations, this effectively increases the statistics 
of the correlator by the number of the grid points as compared to that of the point source. 
Since we compute the quark propagator by splitting it into the low-mode and the high- 
mode pieces 

G = G^ + G^ 

= G^ + $^^.Gf, (30) 

i 

the nucleon correlation function can be split into contributions from the low modes and the 
high modes 



C(G, G, G) = C(G'^ + J2 ^^^f ' G"" + J2 ^jGf, G^ + J^ O^G^) 

i j k 

= C(G^, G^, G^) + J2 G{Q^^\. e,Gt B,G\) 

i 

+ Y, C{9iGf, G^, G^) + J2 GiG"", e,Gf, G^) + J] GiG"", G^, O.Gf) 

i i i 

+ J2 G{OiGt 9,Gt G^) + J2 G{^^Gt G^, e,Gf) + J^ GiG"", e,Gf, 9,Gf) 

i i i 

+ 2_^G{9iGi ,9jGj ,G ) + 2_^G{9iGi ,G ,9jGj) + 2_^G{G ,9iGi ,9jGj) 

i¥=j ij^j i+3 

+ Y. G{9,G\,9,G],9uGl), (31) 

i^i or j^k or k^i 
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where both G^ and ^^ OiGf are noise-estimated propagators. We denote the latter exphcitly 
because we will modify it below to improve the signal-to-noise ratio. 

We first note that the last four terms with summation over different Z3 noises, i.e. terms 
with YliiyLj s-iid Ylii^jorj^kork^i ^^^ pure noise. They can be dropped without changing the 
expectation value of the noise estimate. For the purely high-mode contribution, as well 
as the mixed terms containing two G^ and a single G^, one needs to compute them with 
the Z3 noise technique. However, for the correlators either involving three G^ or with two 
G^ and one G^, one can implement the low-mode substitution technique (LMS) 26| to 
improve the signal. For the case where all the propagators are purely from the low modes, 
we substitute the noise-estimated one (the second term on the right side of Eq. fl3T]) ) with 
the exact expression without noise, i.e. ^,. G{G^ , Gf, Gf ) which is the sum of the low-mode 
contribution to point-source correlators over the grid points. For those terms containing 
only one G^ and two G^, we replace the noise-estimated correlator (terms 6 through 8 on 
the right side of Eq. ( 13T|) with the one where the sources of the two G^ are placed at the 
same grid site i and multiplied with 6f. This ties them with G^ , which contains a 6, to form 
the less-noisy multiple point nucleon correlator over the grid points. (See the upper part of 
Fig. [H) Since we replaced two noise-estimated G^ with one noise-estimated G^ G^, this 
should yield a gain in signal-to-noise ratio by a/Vs, where V3 is the number of grid points. 

Therefore, the nucleon correlator with LMS is 

C^ LMS\G , Gr, Grj = G (Gr , Gr , Gr J + G (G , G , G j + G(G , G , G ) -\- G(G , G , G J 
+ (^ 9^ [G(Gf , Gf , G^) + CiGt G^ Gf ) + CiG'', Gf , Gf )] ) 

i 

+ J2C{Gf,Gf,G^). (32) 



We should point out that this is different from low-mode averaging 3ll, l32| in that the 
noise-estimated low-mode propagators in LMS are substituted with the exact ones whenever 
possible and matched with the high-mode propagators at the grid points, whereas low-mode 
averaging would replace the noise-estimated low-mode part with the average over all space- 
time points. 

Ideally, we would compute a separate set of quark propagators from each source point on 
the grid source, eliminating those gauge-dependent contributions to the two-point function 
(which contribute only noise) where the propagators come from different source points while 
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still gaining the benefits of the extra sources. However, this is prohibitively expensive, as it 
requires additional inversions. As mentioned previously, the grid source is chosen in order to 
achieve some of the benefits of using multiple spatial sources while requiring no additional 
inversions, but noise is introduced by those contributions which mix the various source 
locations. However, the purely low-mode part requires no inversions to calculate, so it can 
be computed exactly without noise from these spurious contributions. Since this low-mode 
contribution dominates the correlation function at large time separation, the use of the low 
modes to compute it exactly allows for a substantial reduction in the error of the nucleon 
mass. 

In Fig. m we plot the effective masses of the nucleon correlation functions from the point 
source, the noise-grid source, the noise-grid source with smearing and folding of the two 
correlators from sources on two time slices at t = and t = 32, and the variation calculation 
with the point and the noise-grid source. The fitted masses are tabulated in Table HI They 
are calculated on the ensemble of 50 2 -|- 1 flavor DWF configurations on the 24^ x 64 lattice 
with ami = 0.005. The overlap propagator is computed with the valence quark mass at 
aTTi ~ 0.016 close to the light sea mass which corresponds to m.^ ~ 330 MeV. We see the 





Point 


Z3-grid 
+ LMS 


Z3-grid + LMS 
+ Smear + Folding 


Z3-grid + LMS 
+ Variation 


Nucleon mass (GeV) 


1.13(14) 


1.08(5) 


1.14(2) 


1.12(1) 



TABLE I. Comparison of nucleon masses from several methods. 

fitted nucleon mass for the point source has a 12% error. The Z3 noise grid source with 
LMS reduces the error by a factor 2.8 down to 4.6%. Replacing the point grid with a 
smeared grid and putting the sources on two time slices reduces the error further to 1.8%. 
Finally, a variational calculation with a second nucleon interpolation field and smeared grid 
sources gives an error of 1.1%. For this, the second nucleon interpolation field we use is 



Xa[X) 



■{u)a/ 



'Xn)b, 



,(d)c 



e "y^^jja ix)tjja {x)C/sj'ijj\ (x) and each quark in the interpolation field can be a 
point-grid or a smeared-grid source. The combination of point-grid and smeared-grid gives 
eight operators. Of these eight, we need to consider the symmetric combination between 
the two u quarks so that they are spatially symmetric. This reduces the choice to six. We 
calculate all of the six combinations for each of two interpolation fields (which differ in their 
Dirac structure) and during the analysis choose four of the twelve which have a good result 
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for the ground state and a decent result for the lowest excited state. The operators of the 
4x4 variational calculation are obtained by choosing x or x' ^"^^ choosing all three quarks 
smeared or just one quark smeared. 

It is worth emphasizing out that replacing the point source by a point-grid source with 
LMS, with no further improvements, reduces the error by a factor of 2.8, improving the 
statistics by a factor of 8 with no added exp ense. For comparison, we note a recent study of 
the all-mode- averaging (AMA) method 33| on the same 24^ x 64 lattices which found that 
to achieve the same error on the nucleon mass with AMA on 32 smeared sources would cost 
1/4 of the computer time of a single smeared source without AMA. This implies that AMA, 
given the same resources, would improve the statistics by a factor of 4. 

The variational approach using two different interpolating fields combined with a mixed 
point/smeared source does reduce the variance, but not enough to justify the increased 
computational cost required for the additional point source inversion. Overall, we find the 
use of smeared Z^ noise-grid sources on two time slices is the most efficient approach; this 
reduces the error of the nucleon mass by a factor of 7 compared to that of the point source 
while adding only the slight overhead of LMS to the computational cost. 



IV. QUARK LOOPS 

The quark loop with zero momentum on a time slice is defined as ^^trrD~^(x, x) where 
the trace is over color and spin and F is a 7 matrix. It is computed for the quark propagator 
which begins and ends at the same site (or more generally a neighboring site if a point- 
split current or smearing function is involved). It is needed for meson correlators with 
the quark-antiquark annihilation channels and the disconnected insertions in nucleon form 
factor calculations. Since the quark loop calculation involves the sum of quark propagators 
originating from all lattice sites on a time slice, it is not practical to calculate it with the usual 
point-source inversion. Instead, stochastic estimation is usually used for such calculations. 
Taking the noise source in Eq. fl20|) where the grid includes all lattice points, the loop at the 
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FIG. 4. The comparison of proton effective masses for the point source, the Z3 point-grid source, 
the Z3 smeared-grid source and variation. The fitted results are given in Table HI 



position i for each gauge configuration is then estimated as 



L(^) = {eiD~\t,x)r^{x)) 






D-H 



(33) 



where the average is taken over the noise. When the noise source r]{x) is diluted in color 
and Dirac indices, the loop L{i) is a matrix in color-spin space. For finite number of noises, 
there is noise contribution from different sites which is not a gauge-invariant quantity and 
is suppressed by gauge averaging in addition to noise averaging. It has been shown that 
Zjv noise is optimal since it gives the minimum variance 3J, l35|. Here we shall use the 
Zi noise which, in addition to calculating loops, can be applied to the evaluation of quark 
annihilation involving two quarks and two antiquarks in the interpolation operators, such 
as in meson-meson scattering calculations. 

In the nucleon three-point function involving the loops in the disconnected insertion 
(DI), the quark loop including the external current needs to be Fourier transformed to give 



16 



a definite momentum transfer. Because of the translational invariance after gauge averaging, 
one does not need to put the noise on all spatial points at a given time slice. Instead, one 
could select an evenly spaced grid separated by A^. sites in each of the spatial directions. In 
this case, the Fourier transformation of f{x) on a lattice with periodic boundary conditions 
can project to a low momentum q in the x-direction via the relation Xljec e~*^^''/(xj). Besides 
the definite momenta q, it involves a mixture of the next higher momentum being qn = 
q±^pi, where L is the spatial lattice dimension and pi is the unit of lattice momentum, i.e. 
pe = ^. Other higher momenta can also be mixed in, but we will not discuss them here. In 
the nucleon DI calculation, the source and sink momenta of the nucleon propagator can have 
a definite momentum (the grid source for the nucleon propagator with A^; = 6 will have a 
zero- momentum source mixed with p = Api). As long as these mixed momenta are taken into 
account in selecting the desired momentum transfer, the mixed qn loop will be suppressed 
due to momentum conservation. Should it happen that the contribution from the mixed high 
momentum loop with qn does not vanish due to the finite noise and gauge configurations, 
it will involve the intermediate state with energy E = \/m1j + qjj. Its contribution is 



exponentially suppressed with a factor ~ e~^^* where AE = \/m?]^] + q^ — a/tti^ + g^. 
Since qn > Q and thus AE > 0, the suppression factor can be substantial for the range of 
t where the signal of the DI calculation is to be obtained. For the present work, the scalar 
content is a forward matrix element which requires g = 0. This makes it easier to consider 
the possible qn contamination. 

In our present calculation of the scalar matrix elements on the 24^ x 64 lattice, we shall 
use the grid with (A^., A^, A^, A^) = (4,4,4,2) for the high modes with odd-even dilution 
as well as dilution in time. This entails calculation of 4 noise propagators (two for odd- 
even and two for time dilution). We show in Fig. [S]a cartoon of the even-odd dilution and 
the time dilution. The nearest neighbor which can give rise to noise contribution is at a 



Euclidean distance of li = v4^ + 2^ ^4.5 which is reasonably far; as such, we do not do 



the unbiased su 
lattice points 



Dtraction as has been done previously for noises which have support on all 
, |36[ , where the nearest neighbor is a single lattice site away. Should there 
be noise contribution from the higher mixed momentum qn despite of suppression due to 
momentum conservation as discussed above, it will be further suppressed by e~^^*. With 
qn = 27r/4, q = and ?Tijv ~ 0.66, AE is 1.04. As a result, the suppression factor for t > 6 
where the strangeness matrix element will be extracted is less than e~^'^^ = 0.002 which 
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FIG. 5. The diagram shows the finite noise contribution of the loops. The red dots are the grid 
sources on even sites. The closed loop shown is an example of the signal and the open-jawed curve 
is an example of a noise contribution which is suppressed by gauge and noise averaging and by the 
finite distance separation of the quark propagator. 



renders this noise contribution negligible. 

To further improve the loop calculation, we use low-mode averaging (LMA) for the low- 
mode contribution by summing over the spatial volume on a time slice, while the high-mode 
contribution is calculated on the grid as described above. 



^ = E^-«+E 



A^AyA, 



LhU), 



(34) 



where Lh is scaled from the sum of the grid points to that of the full space volume to match 
with Liii) from LMA. 

It is interesting to find out the individual contributions to the quark loops from the low 
modes and the high modes. To this end, we measure the following quantity to gauge the 
magnitude of the low-mode contribution relative to the sum of those of the low modes and 
the high modes, 

I T low I 

The norms of the low- and the high-mode parts of the loops are used to avoid the sin- 
gular situation when the full loop is near zero, especially for the vector and axial vector 
channels, where the vacuum expectation values of the loops are zero. The ratio is taken 
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FIG. 6. The plot shows the relative low- mode magnitude vm defined in Eq. (I35|) as a function of 
the accumulated pairs of eigenmodes plus the zero modes. The corresponding pion mass is about 
330 MeV for both the quark in the loop and the light sea quark. The error bars denote the standard 
deviation from 5 gauge configurations. 



before averaging over the time slices so as to be more relevant to the DI calculation where 
the loop is to be correlated with the nucleon propagator in a limited time range. We used 
5 configurations in this study. 

The relative magnitudes of the low modes, defined as tm in Eq. f p5|) . in terms of the 
accumulated pairs of low eigenmodes in addition to the zero modes are shown in Fig. [6] 
for different currents. This is for the case of a quark mass which corresponds to a pion 
mass of ~ 330 MeV. We see that rjv/ of the pseudoscalar loop is dominated by a few lowest 
eigenmodes, and all the higher modes contribute only a few percent. This is not surprising, 
since the space-time integral of the pseudoscalar density times the quark mass is just the 
topological charge of the configuration which is totally determined by the zero modes. For the 
scalar loop, both the L'°" and L'^^^h g^j~g positive on all time slices and for all configurations. 
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FIG. 7. The plot shows the relative low- mode magnitude rM with different bare quark masses for 
the loop. The error bars denote the standard deviation over gauge configurations. This study is 
carried out on 100 configurations. 



Therefore, the relative magnitude tm is the relative low-mode contribution. In this case, 
~ 40% of the scalar loop is contributed by the lowest 50 pairs of eigenmodes. Even at this 
low quark mass, most of the contribution comes from the high modes. In the vector and 
axial vector channels, the relative low-mode magnitude is small with a large variance. 

To study the quark mass dependence, we show rj^ in Fig. [7] with 200 pairs of the lowest 
eigenmodes plus zero modes as a function of the quark mass. The relative magnitude of the 
low modes in the pseudoscalar channel dominates in the light quark region and decreases 
with quark mass. The vm for the scalar, vector, and axial loops decrease with quark mass 
faster than that of the pseudoscalar loop and become very small for m > 500 MeV. This 
hints that for the charmness in the nucleon, the high modes may play a substantial role. 
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V. DISCONNECTED INSERTION OF THREE POINT CORRELATION FUNC- 
TION 

In order to extract the strangeness and charmness contributions to the nucleon given 
in Eq. (E]), we compute the ratio of the disconnected three-point function to the two-point 
function, defined as 

r.u y + ^ < ^^(^) ^^(^') XN{to) > - < XN{t)xN{to) >< qq{t') > /„„x 

""''■ ' ' '°> = < XNm.w > ■ <"" 

The matrix element for the quark content in the nucleon can be extracted from this ratio, 

777 

It, = -^ {N\qq\N) 

777 

= ^ liin R{t,t',to). (37) 

t—t'—>- oo 

For the disconnected insertion, one can write the above ratio in terms of the nucleon prop- 
agator and the quark loop as in Sec. [IV] 



Cit,to)iLit')-{Lit'))) 

R{t, t', to) = ^ ^ '-. (38) 

C{t,to) 



This approach is illustrated in Fig. |8l 

We show in Fig. [9] R{t, t', to) as a function of t' from to = to t = 12 for the case where 
the quark mass in the loop is 0.063, which corresponds to the strange quark mass, and the 
quark mass in the nucleon propagator is 0.016, which corresponds to ?7i^ = 330 MeV. This 
calculation is done on 176 configurations. One can see that the error bars are quite large, 
and it appears to be difficult to determine a plateau region from which to extract the matrix 
element. 



One can improve the statistics by using the summed ratio method 37|]. As illustrated in 
Fig. [ini R{t, f, to) is summed over t' between to + 1 and t — 1 inclusive 

R'{t,to)= Y, Rit,t',to), (39) 

t'=to + l 



Previous calculations of the quark momentum fraction (x) in the DI J38|, have studied 
various choices for the domain of t'. They found that choosing t' between to + 1 and t — 1 
inclusive produced less noise than other choices. For this reason, we have chosen to adopt 
the same domain in our calculations. 
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FIG. 8. The sketch ihustrates the disconnected insertion of the 3pt function. 
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FIG. 9. The plot shows the behavior of R{t,t',tQ) as a function of t' for fixed tg and t at and 
12. The expected plateau in the middle of the time window is rather noisy. 

The ratio method is a means to incorporate information from multiple values of t — to and 
t', thus reducing the statistical error, without needing to explicitly consider which points 
t' might suffer from excited state pollution because they are too close to to and t. The 
method relies on the fact that when the value of R is summed from source to sink, the 
contribution from these contaminated points near the end does not depend on the distance 
from source to sink so long as that distance is large enough. When the propagator length 
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is increased, the additional contribution to the sum conies from a single point in the center, 
within the plateau region. Thus, by examining how much the summed ratio increases when 
the propagator length is increased, we can get an estimate for the plateau value of R that 
incorporates information from multiple values of t' and t. 

that the ratio R'{t,to) has a linear behavior in a region where t is 



It is shown 



37 



large enough to have a plateau in R(t, t', to) and thus 



R'{t,to)^^^mnst.+t{N\qq\N). 



(40) 



We fit the value of the summed ratio R' to constant-plus-linear over the range in which 
t is large enough for R to have a plateau; the slope in R' is thus a measurement of the 
plateau value of R, i.e. the matrix element, from which quark content /r, is determined 
using Eq. (137|) . We should mention that the excited state contamination in the sum method 
is found to be 0(e~^^*) 38|, l39[ where A^; is the energy gap between the nucleon ground 
state and the lowest excited state and t is the sink time of the nucleon propagator (we have 
taken the source time at to = 0). On the other hand, the contamination of the excited state 
in the plateau method is 0{e~^'^^^ + (e~^^(*^~*p))) where tp is the time when the plateau 
appears and ts is the fixed sink time of the nucleon propagator. Since in practice, t ^ tg > tp 
and t ~ ts > (ts—tp), the excited state contamination is less than that of the plateau method 
and this has been demonstrated in the calculation of the isovector qa 39|. 



; i 






(a) (b) 

FIG. 10. The cartoon shows the summed ratio method for the disconnected 3pt function. The 

red part is the additional contribution to the sum when the propagator length is increased, which 

is within the plateau region. 
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VI. NUMERICAL RESULTS 

This work is carried out using valence overlap fermions on 176 configurations with 2 + 1 
fiavors of sea domain wall fermions. The lattice size is 24^ x 64 with the inverse lattice spacing 
1/a = 1.73(3) GeV. The u/d quark mass in the DWF sea is 0.005 which corresponds to a 
pion mass of 331 MeV. The sea strange quark mass is 0.04. This pion mass is matched for 
the overlap fermion with a valence u/d quark mass of 0.016. With the help of the multi-mass 



algorithm 40|, we compute a set of 26 quark masses for the nucleon and 10 for the loops 



with an overhead of about 8' 



'0 



3 



of the cost of inverting the lowest quark mass. The 



summed ratio of the strangeness is shown in Fig. 11(a) for amg = 0.063 and arriud = 0.016. 

The nucleon propagator in the DI of the three-point function is calculated with a smeared 
grid source at t = and t = 32. The time-forward propagator with positive-parity projection 
and the time-backward propagator with negative-parity projection are averaged. We see 
that a linear slope develops after t = 6 where the nucleon starts to emerge in the nucleon 
correlator and it is in the middle of the plateau from Fig. [91 

To see the contributions to R'{t, to) and the respective slopes from the low modes and the 



high modes in the loop, we plot them separately in Fig. 11(b) It is interesting to observe 
that practically all the contributions to -R'(t, to) are coming from the low modes despite the 
fact that the low modes saturate only ~ 15% of the strange quark loop. The high-mode 
contribution is quite small. Since the low-mode contribution is exact, it only has variance 
from the gauge ensemble. The variance of the high-mode contribution comes from both the 
noise and gauge ensembles. Since its contribution is small, there is no need to improve this 
part of the estimate with more noise vectors. 
The renormalized strangeness content reads 

{N\ss\N)n = zf^i2GeV) {N\ss\N), (41) 

where the non-perturbative renormalization constant Z^^{2GeV) = 1.153(5) (stat) has 



been calculated in the RI/MOM scheme 42| with input of Za from the chiral Ward identity. 
It is plotted in Fig. [12] as a function of the valence u/d quark mass. The results are obtained 
from fitting the range from t = 7 to 14 in each quark mass case. The lowest quark mass 
gives ttItt = 250 MeV. 

To obtain a value at the physical point we must perform a chiral extrapolation. The 



mixed-action formula can be derived from the partially quenched expression in 43|, using 
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FIG. 11. (a) shows the plot for the summed ratio of the strangeness. The corresponding pion mass 
is 330 MeV. (b) shows the separate contributions from the low modes and the high modes. 



the simple conversion of partially quenched to mixed-action formulae described in 4^ |45 |. 
The NLO mixed-action extrapolation formula is 46 1 



LO 



{N\ss\N) = {N\ss\Ny'' + 



(4vr^) 



((iV|ss|iV)''^ - (A|ss|A)^^) (J(m^, A,/i) + J(m,„, A,/i)) 



ml 



2 1 

VTiA., — VTiZ 



a'^A„ 



^^•'^> (i^ ^ ^-"'^'ifkF " ^°''" (4./J^ 



(42) 



where (A^, A|ss|A^A) are the leading order contributions to the strange matrix element in 
the chiral limit and the non-analytic chiral loop function is [47| 



J(m, A,/i) = 2A VA2 - m2 -h ze log 



A- VA2 -m'^ + ie\ 



A + VA2 -m2 + iej 



+m log 



m 



/i^ 



-2A2 log 



4A^ 



m^ 



(43) 



where A is the mass difference of A baryon and nucleon. 

In the SU{2) case without an explicit delta degree of freedom, Eq. ( H2l) can be simplified 
to 



ml 



The only mixed action low-energy constant ^mix [44 between the valence overlap fermion 
and sea domain-wall fermions is found to be small J48| - it shifts the pion mass at 300 MeV 
by a mere 16 MeV. We can safely neglect the mixed action effects since they are quite small 
compared to the statistical error. 



25 



Chiral behavior of nucleon strangeness 
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FIG. 12. The dependence of the renormahzed strangeness content on the renormahzed u/d quark 
mass in the nucleon propagator. The round dots show the data points with arus = 0.063 for 
reference and the error band shows the global fitting results with physical strange mass. 



Combining Eq. (HOj) and Eq. (jH]), we obtain a global fit using different time slices and 
quark masses. The fitting model for strangeness reads 



B!{t, to) - [(iV|ss|iV) + A{mi - m°) + 5(m, - m°)] t + Cm„r 



't:^to 



(45) 



where m^ and m^ are the physical quark masses corresponding to the correct vr and K 
masses, and Cmi,ms is a set of constants. 

After the chiral extrapolation in the valence quark mass and the interpolation in the loop 
quark mass, we get the renormahzed strangeness matrix element {N\ss\N)fj = 0.351(65) 
and ms{N\ss\N) = 33.3(6.2) MeV with fitting range from t = 7 to 14 and am'^ = 0.063. 
We also calculate the nucleon mass and extrapolate it to the chiral limit with mN{inT^) = 
'^Af(O) + Ciml + C2ml as shown in Fig. [T^and we obtain m]y{0) = 0.998(39) GeV at the 
physical pion mass. Using this number, we obtain fx^ = 0.0334(62). We plot the recent 
results of ms{N\ss\N) in Fig. [T^ from calculations with dynamical fermions with Nj = 2 
and Nf = 2 + 1 (ETM calculation is with Nj = 2 + 1 + 1). We see that our result has a 



small statistical error and its 5a relative error is comparable to those of Ref. 



M, 



15|, 



16|. 
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FIG. 13. The chiral extrapolation of nucleon mass. 



Similarly, we can compute the charmness content mc{N\cc\N) using the same method. 



The summed ratio of the charmness is shown in Fig. 15(a) and the low/high separation is 



given in Fig. 15(b) We see that the low-mode part still plays an important role, but in 
contrast to the strangeness case, the high-mode contribution is no longer small. 

The chiral extrapolation is shown on Fig. [161 We take am^ = 0.67 which corresponds to 
the physical charm mass in the global analysis of charmonium spectrum with three sea masses 
and two lattice spacings 49|] We find the renormalized charmness content {N\cc\N)fi = 
0.074(24) and /r^ = 0.094(31) after the chiral extrapolation and charm quark mass interpo- 
lation. This number is consistent with the result from MILC collaboration at {N\cc\N) 
0.056(27) 



16|, but with a smaller relative error at 33%. It is interesting to note that 



mc{N\cc\N) at 94(31) MeV is larger than ms{N\ss\N) = 33.3(6.2) MeV and agrees with 



the prediction of 70 MeV in Eq. ( ITOl) from the heavy-quark expansion 

By varying the quark mass in the loop with fixed light quark mass in the nucleon, we can 
check the quark mass dependence of /r^ and /tq. Fig. [T7] and Fig. [18] display the quark mass 
dependence of {N\qq\N)ji and fx^Q respectively as a function of the quark mass. We see 
that the matrix element seems to go down as l/niq at reasonably large m^ and fx^ appears 
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FIG. 14. A comparison of the result of our calculation of the strangeness a term as = ms{\ss\N) 
with those of other groups. The statistic errors are denoted by blue error bars and the total errors 
are denoted by red error bars. 
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FIG. 15. The summed ratio of the charmness is plotted on the left panel. The separate high-mode 
and low-mode parts are plotted on the right panel. These results are computed with a valence 
light quark mass which corresponds to wlt^ = 330 MeV. 
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Chiral behavior of nucleon charmness 
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FIG. 16. The dependence of the charmness content on the u/d quark mass in the nucleon 
propagator. The round dots show the data points with arric = 0.67 for reference and the error 
band shows the global fitting results with physical charm mass. 

to be flat beyond mq ~ SOOMeV. This behavior will be checked with higher precision on the 
32^ X 64 lattice with a"^ = 2.35 GeV which can accommodate heavier quark masses than 
on the present 24'^ x 64 lattice. 

VII. CONCLUSION 

In this work, we have computed both the strangeness and the charmness content of the 
nucleon with overlap valence fermions on 2 + 1-flavor dynamical DWF configurations on the 
24^ X 64 lattice with a~^ = 1.73 GeV and a sea pion mass of 331 MeV. We have employed a 
smeared Z^ noise-grid source with low-mode substitution to calculate the nucleon two-point 
functions, which reduces the error on the nucleon mass by a factor of 7 compared to the 
calculation with a point source. For the loop part of the three-point disconnected insertion 
calculation, we used low-mode averaging to compute the low-mode part exactly and used 
the ^4 noise on a grid to estimate the high mode part. It turns out that the high-mode 
part contributes negligibly to the strangeness content. With the highly-improved nucleon 
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loop quark mass dependence of <N|qq|N>R 
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FIG. 17. The dependence of {N\qq\N)ii on the loop quark mass, with the mass of the nucleon 
valence quarks fixed at a value corresponding to m^ = 330 MeV. 
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FIG. 18. The same as the last figure for the dependence of fx q of the nucleon on the quark 
mass. 
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propagators and the quark loops, we extrapolate to the physical pion mass and obtain the 
precise value fr^ = 0.0334(62) with a better than 5 a signal. The statistical error is quite 
small compared to those of other lattice calculations. The renormalized matrix element is 
(A^|ss|A^)/{ = 0.351(65) and the strange quark sigma term o"s is ms{N\ss\N) = 33.3(6.2) 
MeV. Similarly, we obtain /^^ = 0.094(31) with a 3 a signal which is the first time such a 
result is obtained for the charm beyond a precision of two sigma. The renormalized matrix 
element is (A^|cc|A^)ij = 0.074(24) and the charm quark sigma term is mc{N\cc\N) = 94(31) 
MeV which is consistent with the prediction based on the heavy quark expansion [5|, l50| . 

Even though our present work has a high precision precision, a more meaningful com- 
parison should be with the number of inversions one has to do in order to achieve the 



same precision. To this end, we shall compare with the calculations by Engelhardt [15| and 
JLQCD [18]. They used the direct DI calculation with DWF and overlap which should have 
comparable inversion time as the overlap fermion we use. In our case, we used 176 configu- 
rations and 48 noise vectors to calculate the high-mode part of the quark loop. In the case 

n 

of Engelhardt [15[, ~ 468 configurations were used each with 1200 noise vectors, giving a 
relative error of 24% at sink time T = 10, comparable to our result. However, this approach 

n 

requires ~ 66 times as many inversions as ours. As for the JLQCD calculation [18|, 288 noise 
vectors were used on 50 configurations and the error is 2.54 times larger than ours. Thus, to 
reach the same error as ours, it would take this approach ~ 11 times as many inversions. We 
attribute the efficiency of our approach to the improvement of both the nucleon propagator 
and the quark loop. 

In the present work we have considered the statistical error only. We will continue this 
work on the 32^ x 64 DSDR lattices with a'^ = 1.37 GeV and m^ = 170 and 250 MeV, 
as well as the finer 32"^ x 64 lattices with a~^ = 2.31 GeV and ?7i^ = 290 MeV, in order to 
extrapolate to the continuum limit and the physical sea pion mass in order to address the 
systematic errors. 
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